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Impulse Response Measurement Using Golay Complementary Sequences 
A technique using Golay complementary sequences is introduced for 
measuring impulse responses. Compared to the swept-sine method, it is 
more robust to additive white noise. Free, open-source software is provided 
in Matlab (or Octave) and Pure Data (PD) for carrying out impulse- 
response measurements using the sound hardware found on a typical 
personal computer. 


Introduction 


h(n) 
s(n) r(n) 


Linear system to be 
measured 


[link] depicts a linear system characterized by an impulse response h(n), 
driven by an input signal s(n), and producing the output signal r(n). The 
system identification problem is to estimate h(n) given known 
input/output signals s(n) and r(n). A practical method for identifying finite 
impulse responses uses Golay complementary sequences to excite the 
linear system as described below. Maximum length sequences may be 
alternatively used for this type of measurement, as they are also noise 
signals consisting of 1's and -1's, but perfectly inverse filtering the 
measurement is more computationally intensive. 


Golay Complementary Sequences 


The length L bilevel sequences a(n) and b(n) are Golay complementary 
sequences if and only if the following condition holds, where © denotes the 
autocorrelation operator [link]: 

Equation: 


a(n)Aa(n) + b(n)=b(n) = 2L6(n) 


and 6(7) is the Kronecker delta function. Many references in the audio 
signal processing literature refer to such sequences as Golay codes; 
however, to avoid confusion with Golay error-correcting codes used in 
digital communication, we call the sequences Golay complementary 
sequences. Recall that ({link]) can also be written using *, the convolution 
operator: 

Equation: 


a(—n)*a(n) + b(—n)*b(n) = 2L6(n) 


Given that az (n) and bz (n) are Golay, it turns out that 

aor (n) =[azr(n) bz (n)] and bez (n) = [az(n) —bz (n)] are also 
Golay. This means that Golay complementary sequences can be constructed 
recursively given seed sequences such as ag (n) = [1 1] and 

b(n) =[1 —1]. See the MATLAB/Octave source code generate_golay.m 
for details. Notice also that the resulting bilevel sequences consist of only 
1's and —1's. This means that the signal contains the maximum possible 
power level given that |s(n)| < 1Vn. This property is helpful for 
minimizing measurement noise. 


Let r, (n) = a(n)*h (n) be the response due to the input a(n), and let 

ry (n) = b(n)*h (n) be the response due to the input b(n). Due to ([link]), 
the impulse response h(n) may be determined as follows: 

Equation: 


h(n) = s-(a(n)8ra (n) +b (n)=rs (n)) 


See golay_response.m for more details. 


Measurement Procedure Using Golay Complementary 
Sequences 


. Generate the Golay complementary sequnces golayA.wav and 
golayB.wav using generate_golay.m. 

. Open the pd patch golay.pd, in pd. 

. Ensure that the patch is not in editing mode, and check the “compute 
audio” box in the main pd window. 

. Adjust the “Output Volume” so that when you click on “Record 
Response to Golay A”, the system under test is behaving linearly (i.e. 
not clipping), but so that the input signal to the sound interface is not 
too noisy. 

. If there is an input volume on the sound interface, adjust it so that the 
levels approximately match those shown in [link] when you click on 
“Record Response to Golay A” and “Record Response to Golay B.” If 
the sound interface has no input volume, then you will need to adjust 
the “Output Volume” accordingly. 


golay.pd from the Impulse Response Measurement Toolbox 
RealSimPLE Project 
Edgar Berdahl 
6/10/07 


tabwrite~ golayRespB 


Write Responses To Disk 


write -wave RespA.wav golayRespA 


ite -wave RespB.wav golayRespB 


lsoundfiler 


oOlayRespA 


pd load-tables 


golay.pd after making a measurement with an appropriate input 
level 


. Once you are satisfied with the results, click the “Write Responses to 
Disk” button. 


7. pd will write the files RespA.wav and RespB.wav to disk. Rename 
these files so that the names match the measurement you just made. 
For instance, you might rename them to hpfRespA.wav and 
hpfRespB.wav if they corresponded to the measurement of a high-pass 
filter. 

8. Run golay_response(‘hpf') in MATLAB or Octave to analyze the 


measured response. Plots will be created, and the file hpflmpResp.wav 
will be written to disk. 


Swept Sine Impulse Response Measurement 


Introduction 


h(n) 
s(n) r(n) 


Linear system to be 
measured 


[link] depicts a linear system characterized by an impulse response h(n), 
driven by an input signal s(n), and producing the output signal r(n). The 
system identification problem is to estimate h(n) given known 
input/output signals s(n) and r(n). A practical method for identifying finite 
impulse responses is the swept-sine measurement technique, described 
below. 


Sine Sweep Measurement Theory 


In some cases, it is desirable to relax the power-maximizing constraint 
|s(n)| = 1Vn in favor of obtaining some other desirable measurement 
system properties. For example, we may care more about the accuracy of 
the measurement at lower frequencies compared to higher frequencies, so 
we would like the excitation signal s(n) to contain more energy at lower 
frequencies [link]. We might also be measuring a mechanical or acoustical 
system in which the motor controlled by s(n) behaves weakly nonlinearly. 
If the nonlinearity is memoryless and is NOT preceded by any filtering, 
then the system to be measured matches the Hammerstein model shown in 
[link]. The goal is to measure h(n), independently of the motor nonlinearity 
f(s). Performing the measurement is complicated by the fact that 
superposition no longer holds. 


MEMORYLESS LINEAR 


NONLINEARITY SYSTEM 
oo 
s(n) f(s(n)) r(n) 


Hammerstein Model 


Mathematically, the Hammerstein system behaves as follows: 
Equation: 


r(n) = (f(s)*h)(n) 


It turns out that we can obtain both of these desirable measurement system 
properties by using a new excitation signal s(n). This signal is a sine wave 
with a frequency that is exponentially increased from w to w2 over T’ 
seconds [link]: 

Equation: 


s(n) =sin [K (em — 1) 


I ] 


wy 


where K = “5 and L = i. The MATLAB/Octave code 
wy 


generate _sinesweeps.m generates the appropriate sine sweep. 


Now we consider how to extract the linearized impulse response from a 
measurement. In essence, we need to inverse filter the measurement by the 
excitation signal. 


To this end, we realize that a useful property of s(n) is that the time delay 
Aty between any sample nog and a later point with instantaneous frequency 
N times larger than the instantaneous frequency at s(m9) is constant: 
Equation: 


This characteristic implies that after inverse-filtering the measured 
response, the signals due to the nonlinear terms in f(s) are located at 
specific places in the final response signal. Consequently, the linear 
contribution to the response, which is proportional to h(n) can be separated 
from the other nonlinear terms. We can thus measure a linear system even if 
it is being driven by a weakly nonlinear motor. 


Because the frequency of s(n) increases exponentially, the system is 
excited for longer periods of time at lower frequencies. This means that the 
inverse filter averages measurements at lower frequencies longer, so this 
measurement technique is better suited to especially low-pass noise sources. 


Sine Sweep Measurement Procedure 


1. Generate the sine sweeps using generate _sinesweeps.m 

2. Open the pd patch sinesweeps.pd in pd. 

3. Ensure that the patch is not in editing mode, and check the “compute 
audio” box in the main pd window. 

4. Adjust the “Output Volume” so that when you click on “Record 
Response To The Sine Sweeps,” the system under test is behaving 
linearly (i.e. not clipping), but so that the input signal to the sound 
interface is not too noisy. 

5. If there is an input volume on the sound interface, adjust it so that the 
levels approximately match those shown in [link] when you click on 
“Record Response To The Sine Sweeps.” If the sound interface has no 
input volume, then you will need to adjust the “Output Volume” 
accordingly. 


sinesweeps.pd from the Impulse Response Measurement Toolbox 
RealSimPLE Project 
Edgar Berdahl 
6/10/07 


Write Response To Disk 


write -wave Resp.wav sinesweepsResp) 


lsoundfiler Ipd load-table 


sinesweepsResp 


sinesweeps.pd after making a measurement with an appropriate 
input level 


. Once you are satisfied with the results, click the “Write Responses to 
Disk” button. 

. pd will write the file Resp.wav to disk. Rename this file so that the 
name matches the measurement you just made. For instance, you 
might rename it to nonlinear2Resp.wav if it corresponded to the 
second time you measured the transfer function of a weakly nonlinear 
system. 

to analyze the measured response. This means that the inverse filter 
will be restricted to a dynamic range of 100 (40 dB), which helps 
avoid exaggerating problems beneath w, and above w», where the 
excitation signal has little energy. 0.4 refers to the length in seconds of 


the linear impulse response term to be extracted. Plots will be 
generated, and the file nonlinear2ImpResp.wav will be written to disk. 


Measuring Impulse Responses Using Golay Complementary Sequences and 
Swept Sinusoids 

Two impulse-response measurement methods are demonstrated, with source 
code provided in Matlab and Pure Data (PD). The measurement technique 
using Golay complementary sequences is particularly robust to additive 
white noise, while the swept sine measurement technique is robust to a 
weakly nonlinear motor exciting the linear system being measured. 


Introduction 


h(n) 
s(n) r(n) 


Linear system to be 
measured 


[link] depicts a linear system characterized by an impulse response h(n), 
driven by an input signal s(n), and producing the output signal r(n). The 
system identification problem is to estimate h(n) given known 
input/output signals s(n) and r(n). This module illustrates the swept-sine 
and the Golay complementary sequence techniques for identifying finite 
impulse responses. The first example measures the impulse response of a 
highpass filter using the Golay method, and the second example measures 
the impulse response of a weakly nonlinear loudspeaker driver using the 
swept-sine method. 


Summary Of Objectives 


e To provide a general framework for characterizing single-input, 
single-output linear systems. 

e To demonstrate how to measure the impulse response of a linear 
system using Golay complementary sequences. 


¢ To explain some of the limitations of making measurements using 
standard sound interfaces. 

e To demonstrate how to find the minimum-phase spectrum 
corresponding to a complex spectrum. 

e To explain how to measure the impulse response of a system even if 
the excitation source motor is weakly nonlinear. This measurement 
technique uses a sine sweep test signal. 


Installation 


1. Install a sound card, sound interface, or other full-duplex data 
acquisition card. 

2. Consider testing the data acquisition card viewing the soundcard set-up 
instructions. 

3. Install either MATLAB or Octave. 

A. Install pd. (Alternatively, you may use other software that is capable of 
recording a system's output for a given input excitation signal. The 
software should also be capable of reading and writing WAV files. For 
example, any multitrack recording software should be fine.) 

5. Download tf_meas.zip and unzip the contents into a conveniently 
located local directory. 


Linear System 


Consider the causal, single-input single-output (SISO) system shown in 
[link]. For simplicity, we will take the system to be linear and discrete-time, 
so that it is characterized by its impulse response h(n) or equivalently its 
transfer function H(z), which is the z transform of h(n). 

Equation: 


h(n) & H(z) 


We will assume that both h(n) and H(z) exist so that we can discuss 
measuring them interchangeably. We will further assume that h(n) has 
finite length so that we can measure the response to an input signal in a 
finite amount of time. The goal of this document is to explain how to excite 


the system with a signal s(n), measure the response r(n), and use s(n) and 
r(n) to determine h(n) (and equivalently H(z)). In particular, it is useful to 
pick a signal s(n) that contains a large amount of energy so that 
measurement noise will not significantly corrupt the measurement results. 


h(n) 
s(n) r(n) 


Linear system 


Limitations Of Sound Interfaces 


Sound cards and sound interfaces are not designed for making transfer 
function measurements. They merely provide a cost-effective solution since 
almost all computers have sound cards. We demonstrate these weaknesses 
given measurements made on a PreSonus Firepod sound interface. The 
output from channel 1 was directly connected to the line input on channel 1, 
and the sampling rate was fg = 44. 1kHz. 


1. Sound interfaces attenuate frequencies near DC and near half of the 
sampling rate fs/2. [link] shows the measured magnitude response of 
the sound interface. Ideally it would be constant for all frequencies. 
However, [link] shows that the -3dB point for the DC blocker lies near 
5Hz. The anti-aliasing filters cause the magnitude to roll off sharply 
such that the high frequency -3dB roll-off point is about 20.9kHz. 


Magnitude [dB] 


10° 10° 10° 
Frequency [Hz] 


Magnitude response of the PreSonus 
FirePOD sound interface when the channel 
1 output is directly connected to the channel 

1 line input 


2. The same sound interface filters cause the phase response 
measurement to be incorrect as well. Although ideally the phase 
response would be 0 radians for all frequencies, [link] reveals that the 


phase measurement is distorted in roughly the same regions where the 
magnitude measurement is distorted. 


Angle [radians] 
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Phase response of the PreSonus FirePOD 
sound interface when the channel 1 output is 
directly connected to the channel 1 line 
input 


3. Even when a computer is programmed to pass an input signal flowing 
into the sound interface input out of the sound interface output as fast 
as possible, there is a delay or latency. This delay is typically on the 
order of tens of milliseconds. [link] shows the impulse response 
measured on the PreSonus card. The latency is approximately 62ms. 
The latency could have been decreased some by adjusting software 
settings in pd . The impulse response also rings noticeably due to the 
anti-aliasing filters. 
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Impulse response of the PreSonus FirePOD sound 
interface when the channel 1 output is directly 
connected to the channel 1 line input 


For reference, we provide the measured impulse response shown in 
[link]. 


Minimum Phase Systems 


One further consequence of the delay is that determining the phase response 
of the measured system is more complicated. The delay is responsible for a 
linear phase term since 6 (n — k) + e -?7F*/Js, If the delay is known (or 
measured), then it may be removed by multiplying the measured spectrum 
by e/27fk/fs, However, if the system being measured is known to be 
minimum phase, then the folded-cepstrum method may be applied to find 
the minimum phase frequency response corresponding to the measured 
frequency response. 


The transfer function measurement toolbox assumes that the system being 
measured is minimum phase. This is a valid assumption in many cases. For 
instance, all strictly positive real transfer functions are minimum phase. 
Dissipative systems are strictly positive real (and therefore minimum phase) 
if the appropriate quantity is measured and the sensor and motor are 
collocated. For example, if s(7) controls a motor exerting a force on a 
dissipative system, and r(7) is the velocity at that same point, the 
corresponding transfer function will be minimum phase. This holds for 
other dual variable pairs such as torque and angular velocity, voltage and 
current, and pressure and fluid flow. 


For systems that are not minimum phase, such as systems involving a 
transmission delay between the input and output quantities, the phase 
plotted by the transfer function measurement toolbox is not the system 
phase response, but rather the minimum phase response corresponding to 
the measured system phase response. 


High Pass Filter Measurement 


The circuit shown in [link] was measured using the Golay complementary 
sequence method to show how the sound interface non-idealities affect a 
measurement. V7 was connected to the output of channel 1 of the 
PreSonus sound interface, and Voy was connected to the line input of 
channel 1 of the interface. 
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Vin | Your 
| ‘rd 
High pass filter 


electrical circuit 


The analog transfer function H(f) can be determined analytically using the 
voltage divider rule: 


Equation: 


H(f) = Vour(f) _ R _ gan fRC 


Vin(f)  R+abg  janfRC+1 


In this case, R = 1k and C' = 0. 47 F, so the -3dB point is about 
faaB = ———- ~ 340Hz. [link] and [link] show that the frequency response 
is accurately measured in the range of about 10Hz to about 9 f5/20. 
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Measured magnitude response of the high pass filter 
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Measured phase response of the high pass filter 


The ringing in the measured impulse response distracts from the more 
subtle characteristics of the ideal high pass filter impulse response. For 
transfer functions that pass large amounts of energy at high frequencies, it 
may be more instructive to inspect the frequency domain measurement 
results. 
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Sine Sweep Measurement of a Weakly Nonlinear Loudspeaker 
Driver 


To exeraggerate the nonlinearity of a loudspeaker, we cut the cone of a 
mishandled driver as shown in [link]. We monitored the sound pressure 
several centimeters in front of the dustcap using an Audio Technica 
AT4049a microphone, which has a flat magnitude response to within 3dB 
from 100Hz to 5kHz. The output from channel 1 of the PreSonus sound 
interface was connected to the speaker via a power amplifier, and the 
microphone was connected to the microphone input of channel 1 on the 
sound interface. The following results are typical of sine sweep 
measurements to show how with a weakly nonlinear motor. 


Mishandled driver with 
additional cuts in the cone 


Inverse filtering the measured response results in [link], which is a plot of 
nonlinear2ImpResp.wav. The linear contribution corresponds to the spike at 


the beginning, while the weakly nonlinear terms are clustered closer to the 
end of the response. 
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Full-length response 


The main linear contribution is cut out and plotted in [link]. The 
measurement was not made in an anechoic chamber, so there is a reflection 
about 15ms after the main impact. 
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Linear impulse response term 


The nonlinear terms are shown magnified in [link]. The lower order 
nonlinear terms toward the right have larger magnitude but overlap less in 
time (see [link]). Note that [link] implies that the overlapping could be 
reduced by increasing the total length 7’ of the sweep excitation signal. 
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Nonlinear response terms 


The magnitude and phase responses corresponding to the linear impulse 
response term from [link] are shown in [link] and [link] in blue. For 
comparison, another sine sweep measurement was made at a lower level so 
that the speaker behaved approximately linearly. Decreasing the level also 
resulted in more noise and even some systematic error, as is evidenced by 
the red curves in [link] and [link]. This comparison demonstrates that 
making measurements at larger levels can reduce the effects of noise, while 
nonlinear motor effects can be overcome with the sine sweep measurement 
technique. 
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Proportional to magnitude response of h(n) 
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